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ABSTRACT OF THESIS 


A MOLECULAR-FIELD APPROXIMATION FOR QUANTUM CRYSTALS 

Ground-state properties of quantum crystals have received con- 
siderable attention from both theorists and experimentalists. The 
theoretical results have varied widely with the Monte Carlo calcula- 
tions being the most successful. The molecular field approximation 
yields ground-state properties which agree closely with the Monte 
Carlo results. This approach evaluates the dynamical behavior of each 
-“ir of molecules in the molecular field of the other N-2 molecules. 

In addition to predicting ground-state properties that agree well with 
experiment, this approach yields interesting data on the relative im- 
portance of interactions of different nearest neighbor pairs. Results 
are presented for b.c.c. He^ and b.c.c. He^ at low and high densities 
(down to 10 cm /mole) • Results are also presented for f.c.c. H 2 over 
a similar density range • 
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I. INTRODUCTION 


There have been a number of theoretical Investigations of the 
ground-state properties of quantum crystals over the last decade. A 
quantum crystal is one in which the zero point excursion of a molecule 
from its equilibrium lattice position is a large fraction of the near- 
est neighbor distance. This is in contrast to most crystals where the 
amplitude of the oscillations about .'thfe equilibrium position is small 
compared to the interatomic spacing at low temperature. This allows 
the dynamics of most crystals to be described classically in terms of 
a collection of harmonic oscillators. 

The large oscillations in quantum crystals, a result of small 
atomic mass and the weakness of the attractive part of the interaction 
potential, invalidate the typical classical approach necessitating the 
use of quantum mechanics for their description. In the treatment of 
quantum crystals, another important consequence of the large zero point 
motion must be properly accounted for. This is the strong correlation 
of the motion of neighboring pairs of molecules. The relative motion 
of neighboring pairs is correlated in such a way as to avoid close en- 
counters which involve the strong short-range repulsive forces. 

Most of the recent published calculations employ a variational 
technique with trial wave functions which contain short-range correla- 
tions. The variational approach relies on the Rayleigh-Ritz varia- 
tional principle given by 
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( 1 ) 


where is the ground state energy, H Is the Hamiltonian of the 

system, and is an approximation to the ground state wave function, 

A complete list of symbols appears in Appendix A. The right- 
hand side of equation (1) represents an upper bound to the ground state 
energy which approaches as approaches with equality when 

U) = I'i'o)* 

A simple form for the N-particle wave function which contains 
short-range correlations was first introduced by Jastrow.^ It is 

N N 


^(^^2^,r2> ♦ • • “ I I j i ~ 


( 2 ) 


where is a single particle wave function centered about the 

— N 

equilibrium lattice position R. and the (N - 1) spher- 

1 J K ^ 

ically symmetric two-body correlation functions* The limiting behavior 

lim 

for the two-body correlation functions is f(r) =0 and 

fCr) = 1. Therefore, the f(r) functions lower the probability of 

r-^ 

finding two molecules close together while having no effect at large 
distances. 

The N-particle Hamiltonian is given by 


N 


H 


hi V' 

2m l_j 


N N 

IC=1 ij^j j=l 


(3) 


. th 


where v(r^j) represents the interatomic potential between the i'"" and 

molecules. Evaluation of the right-hand side of equation (1) 

2 

therefore involves the calculation of on the order of N integrals 
all containing 3N dimansions . I'ne differences in the recent theories 
are a result of the manner in which each one evaluates this expression. 
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The first technique used with much success was the work of 

Nosanow^ which employed the use of a cluster expansion similar to those 

used in statistical mechanics. Similar cluster expansion techniques 

3-10 

have been employed in much of the work on quantum crystals. A 

slightly different expansion has been used by Guyer. He and 
Werthamer have each written review articles which cover most of the 
material mentioned above. ’ More recently, Horner has formulated 
a perturbation theory of quantum crystals. 

Most of these formulations give only fair agreement with experi- 
mental data. The cluster expansions, however, converge rapidly enough 
to be useful only for a very limited set of two-body correlation func- 
tions. They require f(r) 'to approach rapidly its limiting value of 
unity. A satisfactory set of correlation functions at low densities 
may be unsatisfactory at high pressures when the nearest neighbor dis- 
tance becomes smaller. An additional problem is that in a cluster ex- 
pansion the completely free variation of both the single particle wave 

functions and the two body correlation functions leads to a liquid 

12 

solution. This is discussed in detail in Guyer 's review article. 

An approach which avoids these difficulties and also comes clos- 
est to matching experimental data is the Monte Carlo representation 

15—19 

used by Hansen and coworkers. Most differences that remain with 

experiment are attributed primarily to an inadequate representation 
for the two-body interaction potential, v(r). 

The work presented here is an important Improvement over the use 
of a cluster expansion. It also removes the restrictions on the two- 
body correlation functicnc; rsce'-^ior;’ for good convergence of the 
cluster expansion. This approach evaluates the dynamical behavior of 
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each pair of molecules in the molecular field of the other N - 2 

molecules. It is found necessary to incorporate in the molecular field 

the approximate motion of the nearest neighbors to the dynamic pair. 

3 4 

This molecular field approximation gives results for He and He which 

agree closely with experiment and the results from the Monte Carlo rep- 

19 

resentation. As in the recent investigation by Hansen the range of 
density covered has been extended to higher densities (10 cm /mole) 
where the cluster expansions used to date are invalid. 

Results are also presented for solid molecular hydrogen. In 
addition to predicting ground-state properties that agree well with ex- 

o ^ 

periment for both H 2 , He and He , this approach yields data on the im- 
portance of interactions of different nearest neighbor pairs, and data 
showing the effects of dynamic motion. 
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II. MOLECULAR FIELD THEORY 

I 

The expectation value for the Hajniltotiian for an N“particle sys- 


tem IS 


<h) ■ (4) 

K=1 j=l 

where the wave function | tj;) is of the form given by equation (2). At 

this point it is convenient to Introduce a procedure first noted by 

20 

Jackson and Feenberg. The expectation value of a typical term in 
the kinetic energy is 


„2 


<♦1 - - 


>1^* 


Using Green's Theorem 


^ 2 — 

j ij) V 1^) dr^^ dr2 • • • dr^j = - 


Vip* • V\ii dr^ dr2 


Therefore equation (5) can be rewritten as 


-ij - V* • v)‘i>^i "^^2 • • ' 


If we allow i|) to be only real then equation (7) becomes 


^K^X^drj_ dr2 . . . dr^ 


by virtue of 


? 1 Z 1 2 

Irrj^ = '7 V ~ --r (Vlp) 


( 9 ) 
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Using a 'f' of the form given by equation (2) gives finally 

N 







^^2 


In t(r .) 
<2 

J 


dr 2^dr^ • • • 


( 10 ) 


We can now rewrite equation (4) as 


<H>- 





" N „ N N 


K=1 


i#J j=l J 


dr^ . . . drj^ 


{j '*’^1 


drj . . . dr^ -1 


( 11 ) 


with 


’13 ^ «C^i 3 > 


(12) 


A common choice for the form of the single particle wave func- 


tion 4>(^jj.) is 


' (f) “P 


“ — (r ^ 
2 


■v=] 


with this form 


(13) 


V^3n4>(r^) = -3g 


and equation (11) further simplifies to 




N N 

i^^j i=l J 


dr^ . . . drj 


(U) 
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and the expectation value per particle is just 


(h> ^ 3^ 

N 4ia 





(15) 


It is at this point that the recent theories for quantum crys- 
tals diverge. Each solves* this equation in a different manner for var 
ious parameterized forms of the two-body correlation function, f(r). 

In the static molecular field theory approximation the effective molec 
ular field acting on an arbitrary pair is found by assuming that the 
N - 2 other molecules are fixed at their respective equilibrium lat- 
tice positions. This is accomplished by letting the Xth and Kth mole- 
cules represent the arbitrary pair and by taking the limit 

^ - R^) for all other N - 2 molecules (i.e., 

i ^ X,tc) . 

Neglecting the constant in equation (15) we can arrange the re- 
mainder so that the sum is of contributions to the energy from suc- 
cessive nearest neighbor shells. Rewriting equation (15) in this 
manner gives 



where the N^'s are the number of neighbors in each shell and (X,k) 

1st n.n. means the Xth and icth molecules are first nearest neighbors, 
etc. dy letting *’• - R.^) in this expression for all 

i 7 ^ X,K we get 



<b) 


4m 


J~'+2(r -r H,2(r^-y£2(r,^) J 7 ^ f^(?- v TT^ ar 


(X,K)lst n.n. 

+ . . . (17) 
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In a previous paper this expression was further simplified by assuming that the products over s 
and p to first order need only be carried out over first nearest neighbors to \ and k. Each func- 

O _ 

tion f (r - R ) was then expanded about its value evaluated at the equilibrium lattice site separation, 
f^(R - R ), These were angle averaged and reduced equation (17) to a collection of easily solvable two- 

K S 

dimensional integrals. 

Recently equation (17) has been evaluated more exactly by calculating the six-dimensional integrals 
numerically using Monte Carlo techniques. The products of two-body correlation functions in equation (17) 
contain all of the correlation functions which couple the Xth and Kth molecules to every other molecule 
in Che system. In practice this product was only extended over the first three nearest neighbor shells 
of Che Xth and Kth molecules. For the two-body correlation function considered in this investigation, 


this was an adequate representation. The values for “ Rg) and 

outside of this range were unity. This approximation is 
easily changed and can be adjusted to include more nearest neighbor 
shells as the two-body correlation functions considered become more 
long ranged. 

Also, the terms in equation (17) extend over the entire lattice. 

In practice only the contributions from the first ten nearest neighbor 
shells are calculated exactly. The contributions from the remaining 
shells are evaluated for a static lattice. This is a more than ade- 
quate approximation. In the Monte Carlo calculations mentioned previ— 
15-19 

ously only the first four nearest neighbor shells are evaluated 

exactly. The details of the evaluation of equation (17) can be found 
in Appendix B. 

The forms of trial wave function and interaction potential used 
in this study are primarily those that have received the most attention 
in other studies. This is because the intent at this point is a compar- 
ison of techniques for evaluating equation (1) rather than a detailed 
study of various trial wave functions and interaction potentials. The 
trial wave functions used in this study therefore are in the form given 
by equation (2), The single particle wave functions in that equation 
are given for this study by' equation (13). The two-body correlation 
function used is the one which corresponds to the WKB solution of the 
two-body Schrodinger equation with a Lennard-Jones potential in the 
limit as r^^ ^0. It has the form 
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Finally the interaction potentials used for this study were of the 
Lennard“^Jones type given by 


v(r^j) = 4ej 

where e = 10.22° K, o = 2.556 S 
a = 2.93 X for hydrogen. 



for helium, and 


e 


37° K, 


(19) 



III. RESULTS AND DISCUSSION 


Equation (17) was evaluated using standard Monte Carlo Integra- 

22 A 

tion techniques for a wide range of variables. The results for He^, 

3 

He , and H 2 were somewhat disappointing. In each case the ground-state 
energy at zero pressure agreed quite well with experiment. However, in 
each case this minimum with respect to our two variational parameters, 

S and <, occurred for a value of 3=0. The corresponding value for 
< was also higher than the value given by investigators using random 
walk Monte Carlo techniques. These results at zero pressure alone were 
not completely surprising since this is near the liquid -solid transi- 
tion. Also, as has been pointed out by other investigators, * the 
ground-state energy is a slowly varying function of the variational 
parameters in the vicinity of the minima. This was the case in this 
study as the value of E^ at 3=0 was not much less than the value 
obtained for the minimizing values of 3 and k reported elsewhere. 
However, the disturbing part of the results was the fact that the mini- 
mum with respect to 3 was at 6=0 for all densities. Even though 
this corresponds to a solid in our calculation unlike the cluster expan- 
sion, it resulted in a trial wave function with effectively only one 
variational parameter, which was undesirable. The resulting ground- 
state energies dropped increasingly below experimental values for in- 
creasing densities. Typical results for He^ are shown in figure 1. 

A. New Formalism 

It -.as felt that these .r»..siijLLS were due to the rigidity of 

the surrounding lattice. The localizing of molecules on individual 
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13 


lattlc. sitaa co.es partisUy fro. the single particle wave functions 
and partially fro. the t.o-hody correlation functions. With surround- 
ing molecules fired on their respective lattice sites, a molecule was 
-onfined to a smaller volume than It would have been If these mole- 
cules were allowed to move. This is another way of saying that the 
correlations with the surrounding lattice were doing more than their 

proper share of the localltlng. yhls In turn allowed the single par- 
tide wave functions to spread out. 

MacMahan calculates the exchange integral for He^, J, using a 
Monte Carlo calculation and also by evaluating the six-dimensional 

obtained by fixins all Kn^ f-Km u 

y rixing all but the exchanging pair of atoms on 

their egulllhrlum lattice sites similar to our approach for the 
g ound state energy. He finds good agreement at zero pressure for 
the two techniques but the results rapidly diverge as the pressure 
increases. He concludes as we do that the rigidity of the lattice is 
to blame and this effect Is more pronounced at higher densities. 

The solution is then to allow In some approximate fashion some 
-vement of the surrounding lattice. It was felt that allowing the 
fctlon of the nearest neighbors of each of the molecules in the dy- 
namic pair (A.r) would be sufficient to remove the difficulties with 
the results. This was verified by later results. 

In an attempt to preserve the dimensionality of the integrals at 
six a number of different expansions were attempted In order to incor- 
porate the approximate motion of these nearest neighbors. Each expan- 
sion was burdened by slow and sometimes questionable convergence. 

Each in turn was abandoned. 

lu older to see clfcmr'v thp i. . 

ft-u-t of allowing some approximate 
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motion to the nearest neighbors of X and k it is convenient to re- 
write a typical integral from equation (17) as 

( 20 ) 

where 

- n U 

sfX,v: pt\,< 

Note that as shown above G(r, ,r ) expresses the correlations of mole- 
cules A and k to all other molecules localized at their equilibrium 
lattice sites. Also expressed in this manner G(r, ,r ) = 1.0 corre- 

A K 

sponds to the cluster expansion result when that expansion is truncated 
after the two-body term. By allowing motion to the nearest neighbors 
to A and k a new G(r,,r ) would be needed which included the ef- 
fects of this motion on: (1) the correlations between the Xth and <th 

molecules and each of their nearest neighbors; (2) the correlations be- 
tween these nearest neighbors with each other; and (3) the correlations 
of these neighbors with all of the remaining molecules localized on 
their equilibrium lattice sites* 

The simplest approximation is to change G(r^,r^) to include only 
the effects on the correlations of the Ath and Kth molecules with their 
nearest neighbors (item 1 above). This gives 
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n 

,K 

?^n.n. to K 


f (r - R ) 

K S 


n 

Pi^X,K 

f n . n • to k 


f^(r. - R ) 

A p 


n 

i=n.n. to X 
3 ^.n. to < 


r - R, 




■(r^ - V^^^'X - . TT / - Rj) 

j=n,n* to K j 
j^n.n. to X 


f^(r 


r . )dr . 
3 ’ 


n S 


^ i=n.n, to 
A and < 


f'(rj - r,)f2(r^ - r^)*2(r^ - R^)dr^ 

( 22 ) 


When molecules X and < do not have any nearest neighbors in common 
this becomes a 48-diTnensional integral for a b-c.c. lattice and a 72- 
dimensional integral for a f^c^c* lattice. Note it is only because we 
have not included the effects of the motion on the correlations between 
these nearest neighbors themselves (item 2 above) that we can write equa* 
tion (22) as the product of independent integrals. 

Evaluation of G(r,,r ) as given in equation (22) was further sim- 
plified by using the same point in three dimensional phase space for the 
evaluation of each of the integrals in that expression. This approxima- 
tion had the effect of adding only three new dimensions to a typical 
integral such as the integral in equation (20) making them nine- 
dimensional integrals . 
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The remaining results are for this new formalism except where the 
original six-dimensional results are discussed for the sake of compari- 
son. This new formalism gave satisfactory results over the whole range 

4 „ 3 

of this study. The following quantities were calculated for He , He , 
and H 2 for the Lennard-Jones potentials mentioned earlier. The ground- 
state energy as a function of density was calculated for each. At each 
density for which the energy was calculated the root-mean-square devia- 
tion from equilibrium was evaluated from 






1/2 


Normal izatlon 


(23) 


Also at each density the single particle distribution function R(r) was 
evaluated. R(t) is given by 

E(|r, - S'J) - "x 

where the normalization R(0) = 1.0 was used. For both 

R(r) , A and < were nearest neighbors. From a curve fit to the energy 

versus density data, pressure and compressibility values were 
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obtained. Some idea of the effect of allowing some movement of the 
nearest neighbors of the dynamic pair can be obtained from figure 2. 

Here R(r) is presented for the two approximations for He^ for two dif- 
ferent densities. Note that the result of allowing some motion to the 
nearest neighbors is a spreading out of the single particle distribu- 
tion function. This is just as expected and is apparently sufficient 
to drastically improve agreement with other ground-state properties of 
solid helium and hydrogen. For the sake of completeness, tabular data 
for the R(r)^s plotted in figure 2 and later figures are all presented 
in the last two tables of the RESULTS AND DISCUSSION section, Tables VI 
and VII. 

B; Helium^ Results 

Results for b.c.c. He are presented in Table I. In this table 

V is the expectation value of the energy per particle from the bare 

Lennard- Jones potential, v(r^^). Also in this table T is the expec- 

2 

tation value of the kinetic energy operator 7^ per particle. For con- 
venience, results for He^ and He^ were obtained only for a b.c.c. lat- 
tice. However, He at high pressures crystallizes in an f.c.c. lattice 
and He^ has the h.c.p. structure for the range of densities studied 
here. This is not very significant since the values of the energy dif- 
fer negligibly for the different lattice structures. Only the values 
of V , and T in the table vary with lattice structure which 

should be kept in mind for any comparisons with experiment or other 
theories - 

Figure 3 contains a comparison between the results of this study, 

19 

py Der imant , and the results of Hansen and Pollack for the ground-state 

2 ^. 

energy. The experimental results are from Pandorf and Edwards down 
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Ground state energy (cal/inole) 








20 


to a volume of 19.5 cm^/mole. Below that, energies were computed by 

integrating the equation of state using pressures from Dugdale and 
26 

Franck. As can be seen from this figure, there is generally good 
agreement with experiment and with the results of Hansen and Pollack. 
This agreement is best at low densities (high specific volume) and the 
calculated energies in this range are within 2 cal/mole of the experi- 
mental values. The statistical errors for the calculated values in 
this range are <1 cal/mole. 

At high densities the calculated values as well as the results 
19 

from Hansen and Pollack both begin to fall below the experimental 
values. This result is also the case for He^ and Is much easier to see 

4 

for He since experimental values exist for the entire density range 

studied. It was concluded in reference 19 that these discrepancies 

with experiment at high density are a result of inadequacies of the 

Lennard-Jones potential which was used, a conclusion which is shared 

15 19 

here. As pointed out in other studies ’ the energy varies slowly 
with respect to the variational parameters in the neighborhood of the 
minimum. Therefore, the values for the variational parameters 8 and 
K should not be taken as exact. Since the values of V and T taken 
separately vary more rapidly with 8 and k these separate values 
have more uncertainty associated with them. It is expected that this 
uncertainty due to not knowing the minimizing B's and <’s exactly is 
less than 10 percent. In comparing the values of V and T with other 
studies keep in mind that these values taken separately will vary with 
lattice structure even though the total energy does not. 

The behavior of the eround energy at high densitie-^ also 

leads to poor agreement with experirueiitai values of pressure and 
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compressibility in this range. Figures 4 and 5 contain comparisons of 
these values. The calculated values were obtained from 


P - 


3E 

— 

3v 


2 

, 3 E 

K = i 2. 

V , 2 

3v 


(25) 


after a simple curve fit of the ground-state energy versus specific 

volume data* The experimental values are from references 26 and 27. 

Figure 6 contains plots of the single particle distribution function 

R(r) for different densities for He . Two of the related values of 

compared with experimental data from reference 24. The 

experimental values are determined from Debye temperature data using 

12 

the Debye formula 


<r^> - I 


i ^ 

4 mK0, 


(26) 


At volumes of 24,5 cm^/mole and 20.8 cm^/mole the calculated values of 
j 3 and 1.06 i, respectively, are both lower than the respec- 
tive experimental values of 1-38 X and 1.15 R. 

4 

C. Helium Results 

Results for b.c.c. He^ are presented in Table II. The format of 
the table is identical to that of Table I. Note that the minimizing 
values of the variational parameters 0 and < are consistently higher 
for He^. A higher B more strongly localizes the individual atoms as 
does a higher value of k. This is physically realistic and is due to 
the heavier mass of the He^ atoms. It leads also to lower values of 

4 

♦■he root- mccn-square : .-i rion end can be seen as a narrowing 

of the single particle distribution functions in figure 7. 


Pressure (atm) 
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Figure 4, Pressure versus volume for He^. 
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Figure 8 contains a comparison of the single particle distribu- 
*8r ^ 4 

tion function and e for He at two specific volumes. The latter, 

_g J.2 

i.e. , e , would be the single particle distribution function if the 
wave function did not contain two-body correlation functions. This 
figure gives some indication then of how much of the localizing of in- 
dividual molecules on lattice sites is due to the effect of the two-body 
correlation functions. Note that at the lower specific volume the two- 

body correlation functions still are responsible for some of the local- 

19 

izing. This supports Hansen and Pollack's contention that these cor- 
relations are still important at high pressures. In addition, this fig- 

15 4 

ure contains a comparison with Hansen and Levesque's R(r) for He at 

a volume of 21.5 cm /mole.* The two R(r)'s appear to be very similar. 

Figure 9 contains a comparison between the results of this study, 

19 

experiment, and the results of Hansen and Pollack for the ground- 
state energy of He^. The experimental data are from references 23 
and 26. Below 11.5 cm^/mole the experimental energies are obtained by 
integrating the experimental equation of state. As for He there is 
good agreement with experiment except at the higher densities (low spe- 
cific volume). As pointed out earlier this is attributed to inadequa- 
cies in the Lennard-Jones potential which was used. An interesting 
feature of figure 3 and figure 9 is that the theoretical values from 
this study and reference 19 both are above the experimental values at 
low densities and are below the experimental values at high densities. 

In each region for both He*^ and He^ the results of this study agree 
more closely with the experimental values. It is not completely clear 
■•'hy this is the case. It Is pcr'':*'’!2 that the use of another lattice 
structure at the high densities would yield a slightly lower ground 






Ground state energy (cal/mole) 
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state energy. This would then lead to closer agreement with the Monte 
Carlo results In the high density region; however, the agreement even 
in this region is remarkably good. 

The differences in the shape of the energy versus density curve 
lead to correspondingly poor agreement with experimental values of 
pressure and compressibility. These results are presented in fig- 
ures 10 and 11 along with experimental values from references 23 and 26. 

3 

As with He , the theoretical values are lower than experimental values 
over the entire range of densities with the worst agreement at the high 
densities. 

The values for the root-mean'^square deviation are not compared 
with experimental estimates" from Debye temperature data since the cal- 
culated values are for a b.c.c. lattice and the experimental values 
would be from a hexagonal close -packed lattice. 

It should be pointed out that even though the Monte Carlo calcu- 
lations of Hansen and coworkers should be fairly exact, there still 

exist differences in the computed ground state energies from Hansen and 
19 IS 

Pollack and Hansen and Levesque. These differences amount to as 
much as 2 cal/mole even though the statistical errors quoted are 
1 cal/mole and 0.4 cal/mole, respectively. 

D; Hydrogen Results 

There exists even more uncertainty in the interatomic potential 
for molecular hydrogen than for the corresponding potential for helium. 
Because of this, the calculations for hydrogen are less extensive than 

for helium. As for helium there exist for hydrogen Monte Carlo calcu- 

23 29 

Iptions well as with which to compare the re- 

sults of this study. The wave fuaccion chosen is the same as the one 
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used for the helium studies. The use of a spherically symmetric two- 
body correlation function f(r) which is independent of the orientation 
of the H 2 molecules requires some justification. Ordinarily one would 
picture the H 2 molecules as something resembling a dumbbell. However, 

30 

calculations by Kolos and Rootaan have shown that the deviation from 
sphericity for the H 2 molecule is only percent. Their calculations 
lead to = 0.84 R and ^ 0.92 X. For a spherical mole- 
cule = ( 22 ^^* Therefore, the use of a spheri- 

cally symmetric two-body correlation function and interatomic potential 
are not as bad as one would at first imagine. 

The results of this study using a Lennard -Jones 6-12 potential 

are presented in Table III,- The ground-state energy at a volume of 
3 

22.65 cm /mole is -170 cal/mole. This Is the same value obtained by 
28 3 

Bruce for a volume of 23-08 cm^/mole and is slightly higher than 
29 

Stewart’s experimental value of -186 cal/mole at a volume of 
22.65 cm^/moie. 

The ground-state energies as a function of volume are compared 
with Bruce’s values in figure 12. The results of this study start to 

3 

drop below Bruce’s values at around a volume of 15 cm /mole and differ 
by over 100 cal/mole at the lowest volumes studied here. The statis- 
tical errors in this region are <10 cal/mole. Most of this difference 
is due to the different parameters used in the Lennard-Jones potential 
in this study. As pointed out previously, these values were e = 37^ K 
and a - 2.93 S. This compares to the values used by Bruce of 
£ = 36.7^ ,K and a = 2.958 R. At these high densities differences in 
the interatomic potential can easily lead to differences of a couple of 
hundred calories per mole in the gLuund-state energy. To show this 
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effect the results from this study have been rescaled to Bruce’s values 
of 0 and £ and are also presented In figure 12. The agreement is 
now excellent at the high densities but not as good in the intermediate 
range . 

Another indication of the amount of uncertainty possible in the 
energies calculated using the Monte Carlo methods of Hansen and co- 
workers is Bruce’s statement that two different initial particle con- 
figurations led to energies that varied by as much as 4 to 6 cal/mole. 
This is in addition to any statistical errors in his results. 

The pressure versus volume results for this study are compared 
with Bruce’s Monte Carlo results and with Stewart’s experimental data 
in figure 13. The pressure, data falls below Bruce’s values at low 
volumes, a direct consequence of the lower energies in this region. 

This brings the pressures from this study closer to Stewart’s experi- 
mental data. This is due to the differences in interatomic potential 
previously pointed out. The resulting pressures are still much higher 
than Stewart's values, a result consistent with all of the theoretical 
studies done to date. Again rescaled values using Bruce's values of a 
and e agree very well with his values for the pressure. Most of the 
studies now in progress on hydrogen are aimed at determining an inter- 
atomic potential which yields theoretical results that fit more closely 
Stewart's pressure data. 

E. Effects of Approximations 

The biggest approximation is, of course, the molecular field 
approximation itself. The effect of studying the dynamic motion of 
pairs of molecules in a lattice '•'’icra most of the remaining molecules 
■ire fixed on their equilibrium lattice sites can only be determined 



Pressure (atm 
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from the comparisons with experiment. This approach in general seems 
to be reasonable based on the agreement with experiment. 

The inclusion of approximate motion of the nearest neighbors to 
the dynamic pair X and < had a dramatic effect on the results. It 
is important therefore to reiterate what this approximation includes 
and what it does not include. This approximation treats exactly the 
correlation of both A and k to each of their nearest neighbors. It 
does not treat exactly the correlation of these nearest neighbors with 
each other or with the remaining molecules in the lattice. This would 
have an indirect effect on the molecules X and k. This effect is 
apparently quite small as Indicated by the close agreement with experi- 
ment with the present formalism. The use of only three new variables 
to evaluate the effects of motion of the nearest neighbors of X and k 
worked quite well again as indicated by the results. 

As pointed out previously, other approximations have been made in 
the evaluation of equation (17). First the integrals were calculated 
exactly only for the first ten nearest neighbor shells. The contribu- 
tions from the remaining shells are evaluated for a static lattice. It 
was found that over the range of densities studied, the energy could be 
determined to within -^l cal/mole if only the contributions from the 
first four nearest neighbor shells were calculated exactly, which is 
the approximation used in the *‘exact” Monte Carlo studies. Evaluating 
contributions from only ten shells exactly leads to a negligible error. 
The static lattice sum is cut off after thirty-eight nearest neighbor 
shells. This again gives negligible error. Some idea of the magnitudes 
of the coritribntions from differunr neighbor shells be seen 

in Table IV, This table contains cuutributions from different groups of 
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nearest neighbor shells for both T and V for He^ at 21-6 cm^/mole, 
and 10*25 cm^/mole. As can be seen from this table, the contributions 
from the first two nearest neighbor shells are very large at the low 
density. The contributions from the other shells become increasingly 
important with higher densities; however, the contributions from the 
first four or five shells dominate the total energy even at these 
higher densities. 

Another important approximation Is the cutoff of the product of 
two-body correlation functions* These products contained the exact 
correlations of the Xth and icth molecules to only their first, second, 
and third nearest neighbors- The values of the correlation function 
outside this range were unity- Some idea of the effect of this approx- 
imation can be seen in Table V where <(t), and are presented 

for a typical calculation at 21.6 cm^/mole as a function of the 

nearest neighbor shells contained exactly in this product- As can be 
seen in this table, these results converge quickly even though the 
f(r) used in this study is a fairly long-ranged two-body correlation 
function. Most of the results for were obtained with a second near- 
est neighbor cutoff in these products co save computation time. Again 
this approximation results in negligible errors. 

The statistical errots ranged from less than 1 cal/mole at the 
low densities to as much as 3 to A cal/mole at 10.25 cm^/mole for the 
helium ground-state energies. The statistical errors for hydrogen were 
a factor of two to three higher. Because of the uncertainties in the 
interaction potential for hydrogen it was not felt warranted to spend 

^ a tii cn t im*^ m o f"* r i > *.5 ^ 
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It is felt that the cumulative effect of errors introduced by the 
above approximations is indeed small and would not effect any of the 
conclusions made from comparisons with experimental data* 

Since most of the results have been compared with Monte Carlo re- 
sults, a discussion of some of the limitations of that approach might 
be useful. These Monte Carlo calculations ^ represented the best 

agreement with experiment to date. However, as pointed out earlier, 

19 15 

calculations by Hansen and Pollack and Hansen and Levesque for 

helium differ by as much as 2 cal/mole even though the statistical 

errors quoted are 1 cal/mole and 0.4 cal/mole, respectively. Also, as 

28 

pointed out previously, Bruce obtained results which differed by as 

much as 6 cal/mole depending on the starting configuration used. This 

was in addition to any statistical error. 

The Monte Carlo method is naturally limited to a finite volume of 

molecules that can be treated exactly. This is usually less than 

1000 particles enclosed in a box with periodic boundary conditions. 

This and other approximations, such as letting each particle interact 

exactly with only four or five nearest neighbor shells, add possible 

errors to the Monte Carlo results, 

31 

In addition, W. W. Wood points out the likelihood that at high 
densities it is quite possible that the Monte Carlo techniques used 
may not properly sample phase space. This might effect substantially 


the results at these densities. 



IV. CONCLUDING REMARKS 


?■ 

¥ 


I 

I 


The molecular field approximation as first formulated (six- 
dimensional integrals) led to undesirable results. The new formula- 
tion which allows in an approximate fashion the motion of nearest 
neighbors of the dynamic pair gives much better results. None of the 
many approximate theories of quantum crystals proposed to date agree 
as well with experimental data and the Monte Carlo results as do the 
results of this study. In addition, limitations on the two-body cor- 
relation functions necessary for the various cluster expansion ap- 
proaches are not imposed by the molecular field approximation. 

The molecular field approximation is, therefore, a valuable tool 
for the study of properties of quantum crystals. The results of this 
study can easily be extended to include different trial wave functions, 
different potentials, and also much higher densities for a thorough in- 
vestigation of the properties of quantum crystals. There are some in- 
dications that random walk Monte Carlo methods may have difficulty in 
properly sampling phase space at very high densities. The molecular 
field approximation might therefore be quite useful at these densities. 
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Volume 

(cm^/mole) 

24.50 

20.80 

16.16 

14.00 

11.82 

10.25 


Volume 

(cm^/mole) 

21.60 

17.50 

15.5 

13.75 

11.82 

10.25 


TABLE I. - b.c.c. He^ RESULTS 


<V> 

(t) 

% 

Pressure 

(cal/mole) 

(cal/mole) 

(cal/mole) 

(atm) 

-44.1 

43.6 

-0.5 


-53.6 

56.4 

2.8 

50 

-64.2 

81.7 

17.5 

260 

-64.0 

102.2 

38,2 

565 

-49.2 

138.5 

89.3 

1380 

-16.5 

183.3 

166.8 

2820 



TABLE II. - b. 

4 

c . c . He 

(v) 

<T> 

Eq 

Pressure 

(cal/mole) 

(cal/mole) 

(cal/mole) 

(atm) 

-54.6 

44.3 

-10.3 


-68.1 

61.3 

-6.8 

93 

-77.3 

77.5 

.2 

250 

-80.4 

97.5 

17.1 

500 

-80.4 

130.8 

50.4 

1035 

-62.1 

165.9 

103.8 

1790 


Compressibility 
(10^ atm-1) 

(S) 

8 

(cj~2) 

fC 

(o"l) 


1.18 

3.5 

1,11 

•23.8 

1,06 

4.1 

1.11 

7.0 

.85 

6.5 

1.08 

3.9 

.75 

8.7 

1.06 

1.8 

.62 

12.5 

1.04 

.7 

.52 

19.0 

1.02 


RESULTS 


Compressibility 

^2)1/2 

6, 

K 

(10^ atm"l) 

«> 

(o 2) 

(o-l) 


1.03 

4.5 

1.13 

9.0 

.83 

5,7 

1.13 

6.0 

.76 

8.0 

1.12 

4.0 

.64 

11.5 

1.11 

2.6 

.56 

15.0 

1.10 

1.6 

.48 

20.0 

1.10 


4N 
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Volume 

(cm^/mole) 

22.65 

18.00 

13.60 

12.00 

10.80 


TABLE III. - f.c.c. H 2 RESULTS 


V 

T 

Eq 

Pressure 

3, 

K 

(cal/iaole) 

(cal /mole) 

(cal/mole) 

(atm) 

(a-2) 

(a-i) 

-301 

131 

-170 


6 

1.19 

-376 

252 

-124 

1 100 

30 

1.14 

-298 

487 

189 

6 350 

67 

1.12 

-51 

595 

544 

12 400 

80 

1.11 

318 

707 

1025 

22 500 

95 

1.10 




TABLE IV. - 

CONTRIBUTIONS FROM 

/ 

DIFFERENT NEAREST NEIGHBOR SHELLS 

Volume 

I^earest 

Contribution to V Contribution 

to T 

Contribution to H 

(cm^/mole) 

neighbor 

(cal /mole) 

(cal/mole) 


(cal/mole) 


shell 

Exact 

Rigid 

Exact 

Rigid 

Exact Rigid 




lattice 


lattice 

lattice 

21.60 

1st 

-29.90 


25.00 


-4.90 


2nd 

-13.27 


4.75 


-8.52 


3rd-10th 

-10.66 


1.94 


-8.72 


llth~38th 

=-0.77 

-0.77 

= .06 

0.06 

=-.71 -0.71 












^V-^=12.51 

4ni 


+12.51 


Total 

-54.60 


44.26 


-10.34 

fS- 

10.25 

1st 

+17.53 


84.90 


+102.43 


2nd 

-38.90 


17.94 


-20,96 


3rd 

-16.37 


3.68 


-12.69 


4 th 

-11 . 03 

-11.38 

2.05 

1.99 

-8.98 -9.39 


5th-10th 

-9.85 

-9.92 

1.40 

1.40 

-8.45 -8.52 


llth-38th 

—3.47 

-3.47 

«.29 

.29 

-3.18 -3.18 





^1^=55.60 

4m 


+55.60 


Total 

-62.09 


165.86 


103.77 



TABLE V. 

- PRODUCT OF 

CORRELATION FUNCTIONS APPROXIMATION 

V 

T 

H 

Exact product of 

(cal /mole) 

(cal/mole) 

(cal/mole) 

f(r)’s cutoff after 

-48.20 

40.80 

-7.46 

1st nearest neighbors 

-50.35 

39.90 

-10.45 

2nd nearest neighbors 

-50.20 

39.90 

-10.30 

3rd nearest neighbors 





TABLE VI 

. - R(r 

) DATA 

FOR He^ 






Volume 






■r(a) 







(cm^/mole) 

0 

0.05 

0.10 

0.15 

0.20 

0.30" 

0.40 

0.50 

0.60 

0.70 

0.80 


10.25 

1.00 

0.928 

0.717 

0.466 

0.254 

0.042 

0.0035 






11.82 

1.00 

.935 

.757 

.559 

.363 

.107 

.019 

0.0018 





14.00 

1.00 

.957 

.824 


.480 

.206 

.060 

.012 

0.0015 




16.16 

1.00 

.961 

.859 


.567 

.293 

.117 

.036 

.0078 




20.8 

1.00 

.975 

.905 


.688 

.445 

.246 

.116 

.046 

0.0148 

0.0039 


24.5 

1.00 


.922 


.738 

.515 

.317 

.173 

.083 

.034 

.012 




TABLE VII. - R(r) BATA FOR He^ 


V olume 






r (o) 





(cm^/mole) 

0.0 

0.05 

0.10 

0.15 

0.20 

0.30 

0.40 

0.50 

0.60 

0.70 

10.25 

1.00 

0.916 

0.681 

0.417 

0.207 

0.027 

0.0012 




11.82 

1.00 

.922 

.730 

.495 

.292 

.064 

.0070 




13,75 

1.00 

.955 

.809 


.412 

.133 

.026 

.0028 



15.50 

1.00 


.829 


.484 

.217 

.068 

.019 

.0022 


17.50 

1.00 


.862 


.587 

.319 

.137 

.045 

.011 

.0019 

21.60 

1.00 


.901 


.675 

.424 

.227 

.102 

.038 

.011 


Six-Dimensional Formalism 


10.25 

21-60 


1.00 0.885 

1.00 .975 


0.610 0.323 0.130 0.0084 0.00015 
.902 .657 .380 .168 0.054 


0.80 


.0028 


0.013 


0.0022 0.00032 



appendix a 


SYMBOLS 


E 

f (r) 
G 
H 
-IT 
i 
J 

j 

m 

N 

P 


R. 

1 


R, . 

^2 


< 




'Jr . 




energy 

two-body correlation function 
defined by eqs. (21) and (22) 

Hamiltonian 

Planck's constant divided by 2 tt 
unit vector in x direction 
exchange integral 
unit vector in y direction 
mass 

number of molecules in the crystal lattice 
pressure 

vector to equilibrium position of the i molecule 
= ^ vector joining the equilibrium positions of the 

i*"^ and molecules 


single particle distribution function defined by eq. (24) 
r distance (depending on usage r = |r^ - r^ | or r = [ r^ - 


etc. ) 


.r2^1/2 


= I r - R I * distance that molecule is from its equi- 

' k: k ’ 

librium position 

root mean square deviation from equilibrium defined by 
eq. (23) 

vector to position of i^^ molecule 
volume ele'Hient for thf’; i loolecule 

= - F. , vector joining the and molecules 
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T kltiGtic energy per molecule 

V potential energy per molecule 

V^j defined by eq. (12) 

V specific volume (appears only in eq. (25)) 

X cartesian coordinate 

y cartesian coordinate 

y = |r^ - R^l and = |F^ % I ’ distance that the nearest neighbors 

to the Xth and xth molecules are from their equilibrium posi- 
tions 

z cartesian coordinate 

z = |r^ - R^|, distance that the Xth molecule is from Its equi- 

librium position 

0 variational parameter in <(>(r) 

6 Dirac 6 function 

V gradient 

2 

V Laplaclan 

E depth of potential well 

6 spherical coordinate 

0^ Debye temperature 

K unit vector in z direction 

< variational parameter in f(r) 

< Boltzmann^ s constant (appears only in eq. (26)) 

K compressibility (appears only in eq. (25)) 

V two-body interaction potential 

7 T angle equal to 180^ 

J.'oLaavje at T.tr.inciLd j. . -i. potential is zero 

$ 


spherical coordinate 
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£j)(r) single particle wave function 

^ total wave function 

it 

complex conjugate of total wave function 
element of solid angle for the molecule 
Subscripts : 

0 ground state 

t h 

1 sometimes associated with the k molecule 

th 

2 sometimes associated with the A molecule 

i 

j 

Jl 


s 

K 

A 

p 


Special symbol: 

expectation value 


APPENDIX E 


DETAILS OF MONTE CARLO INTEGRATION 

In evaluating equation (17) it is first necessary to calculate 

the magnitudes of the following vectors r" - R^.r” - ^ r" - r" 

K k’ X A k’ 

~ ~ sites s and p necessary in 

I. I ^ ~ ^s^ ^ (^»<) pairs evalu- 

sfX,K P^A,K ^ 

ated in the summation. This is done for both b.c.c. and f.c.c. lat- 
tices . 


First we will introduce two new vectors r and z as follows 


r “ r - R 

< K 




Substituting in the other vectors gives 

- R. = F + R 


Ak 


KS 


r, - R = z + R, 

A p Ap 


r -r^=z-r-R, 
A K kA 


(Bl) 

(B2) 

(B3) 

(B4) 

(B5) 


Further let us use the notation that R^ = |r^^| when (A.k) are nearest 
neighbors, = |r^^| when (a,k) are second nearest neighbors, etc. 

First writing r and z in Cartesian coordinates gives 
z = z sin cos + z sin 0^^ sin + z cos 0j^tc, z = [F| (B6) 

r = r sin 0^ cos 4>2^ ^ sin ©2 sin (J) 2 j + r cos ^ 2 '^ » “ |F| (B7) 


With these we can get the magnitudes of the vectors in terms of the six 
variables r, z, 9^^, © 2 , ^ 2 ^ The magnitudes of typical T - R 

K S 

and r^ - used in this study are shown in Table BI for a b.c.c. 
lattice and Table BII for a f.c.c. lattice. Table BIT! cont«in« typi- 
^ lattrces. Tor the later formulation which 
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involved allowing the nearest neighbors of X and k to also move we 
have 

vEr -K = t -R 

s s p p ' 

The vectors r r and r. -- r for this case are similar to those 

K S A p 

in Table BlIIo 

The calculation of equation (17) for a set of variational param- 
eters 6 and k involved evaluating the necessary integrals using 

22 . . 

standard Monte Carlo techniques « This first step involved changing 
the variables of integration so that they ranged from 0^1. Since 
the integrations over r> y, and z ranged from 0 ^ ^ these integra- 
tions were either cut off at a value 10 when step function biasing 
was used or left as integrals between 0 « when some form of expo- 

nential biasing was used. Biasing of the angular variables had very 
little effect and therefore was not used for these variables. 

The procedure used was to pick a point in six-dimensional (or 
nine-dimensional) phase space and then use this same point for a trial 
evaluation of each of the necessary integrals. The 

(- k were evaluated separately. The contribution from 

each nearest neighbor shell involved three separate integrals: one for 

2 

<v(r )>, another for ^ ^ij finally a normalization 

integral. Usually the integrals for the third through the tenth near- 
est neighbor shells were evaluated together since the normalization in- 
tegral was nearly identical for each. Since the contributions from 
these shells are also quite small the effect of this approximation on 
the potential energy was always less than 0.5 percent and on the kinetic 


jLess than 0.1 pL:t:cunt. 
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An exit was made to a statistics subroutine after each set of N 
points in phase space. The program continued until either a set of 
criteria on statistics were met or a given number of points in phase 
space had been evaluated. Typically N ranged from 500 to 1000 and 
the total number of points from 20 000 to 300 000. 
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TABLE BI. - MAGNITUDES OF r - R AND r, - R FOR A B.C.C. 

K S X p 

|7 - R I 

' K s ' 


2 2 ^^^0 

[r + Rq + (sin 02 cos (^2 “*■ ^2 "^2 ^2^^ 

/3 


1/2 


2 , „2 . ^’^^O 


[r + Rq + ;;;j- (-Slll 62 ®2 ^2 '*’ ^2^ ^ 


1/2 


[r^ + R^ + 2rRj^ sin 62 cos 


2 2 1/2 
[r + + 2rRj^ sin ©2 sin <^ 2 ^ 


[r^ + R^ + t/T rR2(sin 02 cos 412 + sin 82 sin 


[r + + /2 rR2(sin 62 <f*2 ^2^' 


1/2 


|r^ - R I 
'X p ' 

2 2 ^^^0 1/2 

[z + R* H (sin 9, cos + sin 9- sin (}i, + cos e,)]"^^ 

U XI i i X 


2 , „2 . ^^^0 


1/2 


[z + Rq + (-sin cos <})^ + sin 9j^ sin + cos 8^^) ] 


2 2 1/2 
[z + Rq + 2zR^ sin 9^^ cos (ji^ J 


2 2 ^/2 
[z + Rj^ + 2zRj^ sin 0, sin 4j^ ] ^ 

[z^ + R 2 + /7 zR 2 (sin 0^^ cos 4>2 "*■ sin 8^^ sin <|)j^) 

(z^ + R 2 + /r zR 2 (sin 6^^ cos + cos 0^) 


LATTICE 
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TABLE BII. - MAGNITUDES OF r - R AND r, - R FOR A F.C.C 

K s X p 

|r - R I 

' K s' 

+ Rq + rRQ(sin 6 ^ cos 4*2 sin ©2 sin 
[r^ + Rq + /J rRQ(sin 8 ^ sin 412 + cos 62)!^^^ 

[r^ + r| + 2rRj^ sin ©2 cos 

[r^ ■*■ 2rR^ sin 02 sin 

It". - R I 
'a pi 

[z^ + Rq + zR^Csin 8 ^ cos + sin 9j^ sin 

[z^ + Rq + zRQ(sin sin + cos 

2 2 1/2 
[z + + 2 zRj^ sin cos ' 

[z^ + R^ + 2zRj^ sin 6^ sin 
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TABLE Bin. - MAGNITUDES OF r., FOR A F.C.C. AND B.C.C. LATTICE 

ak 


b .c^c« , r 


Ak' 


.\e.irest neighbor 
pair 


2 . „2 . 


2zRr 


1 st; [z + r + Rq - 2 zr(sin 6 ^ sin 0 ^ cos((frj^ - (^ 2 ) + cos 6 ^^ cos 62 ) (sin 0^^ cos 


2rR. 


. . '0 
+ sin 0j^ sin + cos 6j^) + ■, (sin 02 cos ^2 + sin 02 sin 412 + cos 02 )] 


1 st 


2 . .2 , „2 


[z + r ■*■ ^ 2 ^ ~ 2 zr (sin 0^^ sin ©2 cos(tj>^ - <^ 2 ) + cos 0 j^ cos 02 ) - 2 zR^ cos 0^ + 2 iR^ cos 62 ] 


1/2 


f .c.c. , r 


Ak ' 


[z^ + + Rq 


+ r2 r2 


2 zr(sin 0 ^^ sin 02 cos(({)^ - 42 ) + cos 0^^ cos 02-/2 zRQ(sin 0^^ cos 

+ sin 0j^ sin + /l rRQ(sin 02 cos <J >2 + sin 02 sin <{> 2 ) ]^^^ 

2zr(sin 0^^ sin ©2 cos((|)^ - <(>2) + cos 0^ cos © 2 ) - 2zR^^ sin 0^^ cos 


Ui 

Lo 


+ 2rR^ sin ©2 cos ] 


1/2 


2nd 


